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ABSTRACT 

The statistical properties of a complete, flux limited sample of 197 long gamma-ray bursts (GRBs) detected 
by BATSE are studied. In order to bring forth their main characteristics, care was taken to define a representa- 
tive set of ten parameters. A multivariate analysis gives that ~ 70% of the total variation in parameter values 
is driven by only three principal components. The variation of the temporal parameters is clearly distinct from 
that of the spectral ones. A close correlation is found between the half-width of the autocorrelation function (t) 
and the emission time (750); most importantly, this correlation is self-similar in the sense that the mean values 
and dispersions of both r and 750 scale with the duration of the burst (Tgo). It is shown that the Amati-relation 
can be derived from the sample and that the scatter around this relation is correlated with the value of r. Hence, 
r has a role similar to that of the break in the afterglow light curve (ft,) in the Ghirlanda-relation. In the standard 
GRB-scenario, the close relation between a global parameter (ft,) and a local one (t) indicates that some of the 
jet-properties do not vary much for different lines of sight. Finally, it is argued that the basic temporal and 
spectral properties are associated with individual pulses, while the overall properties of a burst is determined 
mainly by the number of pulses. 

Subject headings: gamma rays: bursts - gamma rays: observations - methods: data analysis 



1. INTRODUCTION 

The Burst and Transient Source Experiment (BATSE) col- 
lected during its operation the most extensive gamma-ray 
burst (GRB) database to the present day, with more than 2700 
triggered bursts. Even with the advent of Swift and the fu- 
ture mission GLAST, it will remain the most complete burst 
catalog for many years to come. With such a wealth of data, 
there are many statistical techniques that could be used to in- 
vestigate the sample characteristics. So far this analysis has 
focused mainly on bivariate methods. Such an approach has 
some drawbacks for intrinsic relations involving more than 
two observables, since the analysis is limited to two dimen- 
sional projections of the actual relation. Not only do such 
projections increase the scatter but, most importantly, selec- 
tion effects and observational biases can, for example, intro- 
duce spurious correlations. Hence, multivariate methods are 
to be preferred when the quality of the data is such that this 
can be reliably done. In the present work a principal com- 
ponent analysis (PCA) is performed on a complete, flux lim- 
ited sample of GRBs chosen from the BATSE database. Al- 
ter native multivariate meth ods have been discussed and used 
by Mukheriee et al. ( 1998) in connection with attempt to dis- 
criminate between distinct classes of GRBs. 

The PCA is a well-known statistical tool, widely used in 
many scientific disciplines. Some areas of astronomy, like the 
study of active galactic nuclei and quasi-st ellar objects, have 
long benefited from this technique (e.g., by Boroson & Green 
fl992l lB oroso nl2002l) . The method finds, in decreasing order, 
orthogonal directions of maximal variation. The most com- 
mon use of the PCA is to reduce the number of variables 
needed to reproduce the variability of the data without es- 
sential loss of information. Since in most practical situations 
there is a significant amount of inter-correlation among the 
observed quantities, there is then some degree of redundancy 
in the multi-dimensional data. This simplifies the data analy- 
sis, for example, when looking for subcl asses. Previous uses 
of the PCA in the field of GRBs (e.g., by Bagoly et al.ll998t 



Balazs et1ifll2"00l were done mainly for this purpose. How- 
ever, the PCA can also be an important aid to better under- 
stand the generally complex correlation structure encountered 
in multivariate studies and to gain physical insight to the sys- 
tem under study. It is with this purpose in mind that we have 
selected a broader range of physical parameters for our sam- 
ple of GRBs than the ones considered in previous works. We 
include various basic temporal and spectral parameters in or- 
der to characterize temporal variability and spectral evolution. 

Throughout this paper only long GRBs will be considered 
(i.e., those with time duration > 2 s) since short bursts need to 
be treated separately. This is not only because they are most 
likely a different physical class, but also because it would not 
be possible to derive for them some of the parameters under 
consideration, which are already hard to measure for the long 
bursts. The characterization of the temporal properties is com- 
plicated by the remarkable morphological diversity of GRBs 
light curves and the presence of long emission gaps in a sig- 
nificant fraction of them. Hence, a parameter designed to give 
a measure of some property for a "typical" GRB light curve 
will often fail to give consistent or physically meaningful re- 
sults in other cases. The problem extends to any other tempo- 
rally averaged property which follows the pulse structure of 
the light curve. For this reason, seemingly similar measures 
of a given property but based on different criteria can have a 
complementary roll, as it will be shown in this work. 

The paper is organized as follows. In §|2]we present the 
sample and discuss the motivation for the choice of parame- 
ters used in the subsequent analysis. A short introduction to 
PCA is given in § [3] Special attention is given to the esti- 
mate of uncertainties in the coefficients of the principal com- 
ponents; in particular, the effects of mixing and reordering of 
the principal components due to (nearly) degenerate eigenval- 
ues are discussed. The basic statistical results are presented in 
§ @] where also a brief discussion is given of how they com- 
pare to previous studies. A more in depth analysis is done in 
§ |3 where it is shown that two redshift independent correla- 
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tions can be derived from the data. One of the challenges of 
statistical studies is to distinguish between direct causal con- 
nections and those due to association through some common 
variable not necessarily included in the analysis. In an at- 
tempt to make the overall correlation structure o f the parame- 
ters more distinct, we develop a new method in § 15.31 in which 
ranking of correlations is combined with the PCA. It is shown 
that this reveals considerably more structure than does, for ex- 
ample, the use of only the ranking of correlations. Finally, in 
§Hjwe discuss our results with a particular emphasis on how 
they can be used together with results already established to 
illuminate some of the basic properties of GRBs. 

2. DATA 

2.1. Burst Data - Sample Selection 

The present work is based on data taken by BATSE 
on board the Comp ton Gamma-Ray Observatory (CGRO; 
Fish maiTet alJ [l989) that operated between the years 1991— 
2000. It consisted of eight modules placed on each corner of 
the satellite, giving full sky coverage. Each module had two 
types of detectors: the Large Area Detector (LAD) and the 
Spectroscopy Detector (SD). The former had a larger collect- 
ing area and is suited for spectral continuum studies, while the 
latter was designed for the search of spectral features (lines). 
The CGRO Science Support Center (GROSSC) supplies a 
number of high-level data products in its public archive. For 
our spectral analysis we used primarily the high energy reso- 
lution (HER) background and burst data types from the LADs 
having 128 energy channels in the 20-2000 keV energy band 
summed from triggered detectors, and a time-resolution of 
multiples of 64 ms. For long and bright bursts it is often 
the case that the recorded HER data are incomplete, however 
for the bright cases we were able to used the Spectroscopy 
High Energy Resolution Burst (SHERB) data obtained with 
the SDs. These detectors had better energy resolution (256 
energy channels) in the same energy band, but given their 
smaller collecting area, the signal-to-noise ratio (S/N) is of- 
ten too low for time resolved spectral analysis. If that was 
the case, we used alternatively the Medium Energy Resolu- 
tion (MER) data that consist of 16-channel spectra. When 
possible, comparisons of the best-fit parameters using these 
three data formats were made, and the results were consistent 
within uncertainties (see also Preece et al. 2000). 

For the temporal analysis we used the so-called concate- 
nated 64 ms burst data, which is a concatenation of the three 
standard BATSE data types DISCLA, PREB, and DISCSC. 
All three data types have four energy channels (approximately 
25-55, 55-1 10, 1 10-320, and > 320 keV). The DISCLA data 
is a continuous stream of 1.024 s and the PREB data covers 
the 2.048 s prior to the trigger time at 64 ms resolution, both 
types obtained from the 8 LADs. 

To construct our sample, we first selected from the Current 
BATSE Catalog all bursts with time duration Tgo > 2 s and a 
peak flux measured within the 1 .024 s time scale Fi s > 4 pho- 
tons cm" 2 s _1 in the 50-300 keV energy band. A few cases 
(11) were not suitable for our analysis due to, e.g., data gaps, 
poor background information, or missing data formats, leav- 
ing a total sample of 197 bursts. 

2.2. Parameter Selection 

In order to characterize the total time of activity of a burst 
different criteria can be adopted. The most commonly used 
in the literature is 790, the time interval within which 5% 



to 95% of the total count fluence has been detected. Simi- 
larly 7go measures the interval during which 50% of the total 
counts have been observed. In general, these two quantities 
are highly correlated, but the second is sometimes preferred 
when studying "weak" bursts as a more robust measure, but 
this difference is not relevant for our sample. We used values 
taken from the BATSE catalog. 

Since a significant fraction of bursts show quiet intervals 
or gaps in their time histories, other complementary parame- 
ters have been proposed to estimate the emi ssion time scale. 
We us e the emission time T„ as defined by Mitr ofanov et alJ 
(1999), which is the smallest time interval that contains a 
given percentage n% of the total count fluence. Again, 50% 
and 90% are the most commonly used percentile values, and 
we choose here 7go since it is less correlated with 790 (i.e., 
it provides less redundant information). To calculate it, first 
the bins of the (background-subtracted) count time history are 
sorted by intensity. Then, starting from the largest intensity 
bin, 750 is the integration time needed to accumulate 50% of 
the total counts. 

We considered various parameters related to the temporal 
variability of GRBs that have been previously used in the lit- 
erature. The parameter used here is the half- width at half- 
maximum of the autocorrelation function (ACF), which will 
be denoted by r, and it was calcu lated in the 50-300 keV 
energy band following Borgonovo (2004). In the context of 
GRBs, where most light curves show multiple uncorrelated 
pulses, r gives a measure of the typical pulse time-scale. 

A dimensionless variability parameter, calculated in the lo- 
cal frame of the burst, was proposed by Reic hart et alJ ll2001l) 
as a possible "Cepheid-like" parameter that would correlate 
with the burst luminosity (though see Guidorzi et a OHM 
for a critical analysis of this claim). Since, in the present 
analysis, we are considering the observed characteristics of 
the bursts (i.e., at redshift z = 0), we will a dopt a some what 
simplified version of that given by Reichart et al. (2001J). For 
a uniformly sampled light curve, let c, be the net count level 
at bin i, and fo, the corresponding background level. The vari- 
ability V is then defined as 

v= (£ [fa- *) 2 -(c- +*/)]) 

where s, is the corresponding value of the smoothed light 
curve using a boxcar function with a width given by the 
emission time at the 45% level, 745. Thus, the parame- 
ter is a measure of the variability content at high frequen- 
cies. Subtracting the gross detected counts in the numera- 
tor takes into account the Poisson noise contribution, while 
the sum in the denominator is used for normalization. Al- 
ternatively, we could have used the square of the number 
of counts at the peak of the light curve, FA, as first pro- 
posed in Fenimore & Ramirez-Ruiz (2000). In a recent study 
iLi & Paczvfekil "( l2006l) found a better variability-luminosity 
correlation using the F p 2 k normalization, although they also in- 
troduced a different filtering method. Initially we considered 
both options, but in our case the latter choice shows similar 
though weaker correlations. Therefore, we will only present 
the results obtained from the variability definition given in 
equation Q. The percentil e level of the smoothing time- 
scale 745 was determined by Reichart et al. (2001) to render 
the best results in their variability-luminosity correlation stud- 
ies, which gives an additional motivation for the inclusion of 
the nearly identical 750 in our analysis. 
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TABLE 1 

Statistic Description of Burst Parameters 



Parameter 


Median 


QD 


Min 


Max 


<r(logA") 


K-S test 


G test 


7k)(s) 


31.2 


22.7 


2.24 


674. 


0.48 


0.42 


0.57 


T 50 (s) 


4.67 


3.63 


0.256 


32.6 


0.44 


0.11 


0.03 


r(s) 


3.12 


2.81 


0.144 


39.2 


0.49 


0.70 


0.29 


V 


0.0039 


0.0030 


0.000084 


0.029 


0.49 


0.018 


0.32 


St 


1.49 


1.33 


0.086 


57.3 


0.52 


0.29 


0.34 


T)ag(s) 


0.053 


0.034 


0.001 


2.70 


0.51 


0.0002 


0.0003 




1.36 


0.40 


0.40 


6.05 


0.17 


0.005 


0.42 


•FflO 6 erg cm" 2 ) 


35.3 


35.0 


1.70 


781. 


0.54 


0.31 


0.18 


£ pt (keV) 


35 1 . 


146. 


70. 


2251. 


0.28 


0.76 


0.11 


a 


-0.59 


0.27 


-1.59 


1.49 




0.001 


< 10" 5 



The emission strength of the bursts is usually characterized 
by either its peak count flux in a given time interval or its en- 
ergy fiuence, both directly available from the BATSE catalog. 
Since the peak flux F\ s has been used for our sample selec- 
tion, its observed distribution is more limited and biased than 
in the fiuence case. Restricted to the sample range, the peak 
flux shows weak or no significant correlation with the other 
parameters considered in our PCA, with the exception of the 
fiuence. For this reason in the presented PCA we will only 
include the total energy fiuence T. Fiuence values discrim- 
inated by energy channel were first considered but they do 
not show significan t differences. In their correlation analysis, 
Bagol vet alJ {l])98) found the fiuence from the fourth energy 
channel (> 300 keV) to be individually a significant principal 
component. Note, however, that their sample comprise both 
short and long bursts. When considering long bursts, only 
the principal component given by the sum of all energy chan- 
nels was found to be significant, i.e., the fourth channel fiu- 
ence discriminates between the two well-known soft/long and 
hard/short classes but does not stand out within each class. 

To characterize the GRBs spectral shape we modeled the 
photo n spectrum with the standard GRB function ( Ba nd et all 
1993), which is a softly broken power-law. This empirical 
function has four free parameters; namely, the amplitude, the 
peak energy £ p k at which the fifE-spectrum is at its maxi- 
mum, and the low and high asymptotic power-law indices a 
and (3. Unfortunately, in general the high-energy index j3 is 
poorly constrained, and often it has to be fixed to allow the 
convergence of the modeling algorithm. Hence, we included 
in our study only the spectral parameters E p ^ and a. Since 
all long GRBs show a significant spectral evolution, initially 
we considered both time integrated spectra over the whole 
burst and time-resolved spectra at the time f max when the light 
curve reaches its maximum. Note that because we required a 
S/N > 40 to allow a reliable spectral deconvolution with the 
detector response, the time resolution at the light curve peak 
varies for different bursts, but it is typically of ~ 128 ms. The 
spectral parameters measured for the two cases are highly cor- 
related. Since the correlations with the other variables were 
strongest for the f max -case, only this case is included in our 
study 

Except for the simple bursts with a single long smooth 
pulse, where the characteristic energy changes monoton- 
ically from hard to soft, most bursts show a very complex 
spec tral evo lution that correlates with the pulse structure (see, 
e.g., iBor^onovo & Rvdel 120011) . To have a simple measure 
of the overall change, we first divided the burst light curve 
into two fluence-halves, from 5-50% and 50-95% of the total 



counts, and then determined the individual integrated spec- 
tra for each of the two time intervals. We defined the ratio 
TvLEpk = £pk (1 '/Epk (2) > so that TZE V ^ > 1 implies an overall hard 
to soft evolution. The difference between the a-values for 
each fluence-half was initially also considered, but in general 
the estimated changes were not very significant when com- 
pared with the large a uncertainties. This results in a very 
noisy parameter that only shows weak correlations with the 
others, and it was therefore excluded. From the division in 
two fluence-halves we can derive another related parameter 
by considering the ratio between the corresponding time in- 
tervals. Thus, we define the emission symmetry as = T2/T1 
so that Sf > 1 implies an overall decay of the intensity on 
the pulses, or in the case of single-pulse bursts, a fast rise and 
slow decay. Note that by definition it follows that Tgo = T\ +T2. 
The two new parameters 1ZE p k and have obviously the im- 
portant property of being, at least to first order, redshift inde- 
pendent. 

A more commonly used parameter related to the spectral 
evolution of the bursts is the time-lag between two energy 
channels, derived from th e analysis of the cross-correlation 
function (CCF). Following Norris (2002), we define the time- 
lag, 7~i ag , as the time at which the CCF reaches its maximum. 
Among all possible CCFs between the different energy chan- 
nels, we consider here only the one between the third a nd the 
first (CCF31) that gives the best S/N jNorris etalT2 000). Note 
that the parameters TZE v y_ and Ti ag are in principle complemen- 
tary, since they are related to the spectral changes on long and 
short time-scales, respectively. Therefore, in general, the first 
is a measure of the overall spectral evolution and the other of 
the evolution within a pulse. 

In summary, from all the GRBs parameters considered in 
this section, we have selected 10 for multivariate analysis, de- 
scribing a broad range of temporal and spectral characteris- 
tics. They are the following: duration time T^, emission time 
750, ACF half-width r, variability V, emission symmetry 
CCF31 time lag ri ag , the ratio of peak energies HE V ^, fiuence 
T, peak energy E p ^, and low energy spectral index a. 

3. METHOD 

For a data set with n variables, the most basic way to an- 
alyze the multidimensional structure of the correlations be- 
tween all the variables is to organize the linear correlation co- 
efficients R between all possible pairs of variables in matrix 
form, i.e., to construct a correlation matrix. However, unless 
the resulting matrix is very simple or the number of variables 
very small, usually not much direct information will be gained 
by this approach. Various methods have therefore been de- 
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veloped with the aim to simplify the problem enough for the 
main characteristics of the data to become discernible. The 
principal component method achieves this essentially by di- 
agonalizing the correlation matrix and consequently provides 
a new set of uncorrected variables (called principal compo- 
nents, PCs) associated with the simplest possible matrix cor- 
relation structure (i.e., the diagonal form). 

Principal component analysis (PCA) is a widely used tech- 
nique in multivariate analysis. By definition, the first PC gives 
the direction of maximal variance. The second PC gives the 
direction of the second highest variance orthog onal to the first 
PC, and so forth. It has been shown (see, e.g., Jolliffe 2002) 
that the computation of the PCs reduces to the problem of 
finding the eigenvectors and eigenvalues of the correlation 
matrix derived from the original variables, where the eigen- 
values give the variance associated with each PC. Real num- 
ber solutions are ensured by the matrix symmetry. The PCs 
are then sorted in decreasing eigenvalue order, so that each 
PC has associated a decreasing percent of the total variability, 
given by the n-dimensional variance of the sample. They are 
usually defined with unity norm so that, except for an arbitrary 
sign, they are uniquely determined. In theory, this is not gen- 
erally true since there could be degenerate eigenvalues, but 
this situation is extremely unlikely when dealing with experi- 
mental data. However, the existence of very similar eigenval- 
ues can have practical consequences, as we will discuss later 
in connection to the uncertainty estimations. 

By construction, every PC is a linear combination of the n 
original variables Xj 

n 

PCi = ^aijXj (2) 

7=1 

where the coefficients a,y are sometimes referred to as load- 
ings since they represent how much each variable Xj con- 
tributes to a given PC,-. By definition, the matrix of coefficients 
A = (ajj)„ x „ transforms the data points expressed in the vec- 
tor space basis comprised by the original variables Xj to the 
PCi basis, and it is obtained by arranging the n eigenvectors 
of the correlation matrix in columns. Consequently, each row 
in A is associated with one PC and each column with one of 
the original variables. 

The PC representation has many advantages, and it is most 
commonly used to reduce the number of variables needed to 
described the data. Since the first PCs retain most of the vari- 
ability of the data, we could then discard the last PCs with- 
out significant lost of information. Various criteria have been 
proposed to decide how many PCs should be reasonably kept. 
Once the number of variables are reduced to a few, the anal- 
ysis of the data and their graphical representation are greatly 
simplified (e.g., when looking for subclasses or clusters, and 
classification). A complementary aspect of the PCA is that 
it can be used to find near-constant linear relationships be- 
tween the variables when the last PCs have small associated 
eigenvalues (and therefore very small associated variances). 
In those cases, the interpretation of the last PCs is conse- 
quently much more direct than that of the first PCs, and it 
can often be of considerable interest and provide constrains 
on physical models. 

We will use the logarithms of the variables defined in § 12.21 
since the distribution of most of our variables are well de- 
scribe d by lo g-no rmal distributions (see, e.g., iPreece et alJ 
2000, and § 14.11 below) and taken into consideration that 
most of the inference based on PCA assumes a multivariate 



normal distribution dJoHiffel2002D . Thus, unless stated other- 
wise, throughout this paper we will speak about linear corre- 
lations but between logarithmic variables. The only variables 
in our set that are not positive by definition are a and T\ ag . For 
the low energy index, we considered both the direct value of 
a and the simple transform log(a+2) since a+2 > for all a, 
but the results were very similar, and only the first option will 
be presented. In the lag case, 3 points in our sample have near- 
zero negative values, but within errors they can be take n as in- 
trinsically positive. For a much larger sample. iNorri si J2002I) 
also found that the lower tail of the lag distribution reaches 
negative values, but he has shown that the observations are 
consistent with the assumption that all GRB lags are positive, 
and the few exceptions are produced by measurement uncer- 
tainties. We considered both the use of logTi ag with the ex- 
clusion of the 3 negative points, or the adoption of positive 
values within uncertainties (la), with almost identical results. 
The last option was chosen for the present analysis. 

To estimate the uncertainties of the PC coefficients we used 
the bootstrap method (see, e.g. JPress et all 19 92) which does 
not rely on any assumptions of the distributions. However, the 
use of this method in the PCA is not straightforward. To ob- 
tain meaningful results, for each generated surrogate sample, 
care must be taken both to prevent sign switching (since they 
are arbitrary) and reordering of the PCs. The last problem 
occurs when two eigenvalues are equal to within uncertain- 
ties. This situation will often produce a switch in order, but, 
in practice, it may act also as a degeneracy where linear com- 
binations of the two PCs could result. For this purpose, our 
code implementation of the bootstrap method compares the 
original matrix of coefficients A with each bootstrap realiza- 
tion B. By construction A is orthonormal, i.e., the product 
with its transposed is AA' = I, where I is the identity matrix. 
Similarly, the product matrix C = AB' has most elements in 
each row nearly zero, and one element (or more in case of 
degeneracy) approximately ± 1 . After identification of the ex- 
treme elements in each C row, the PCs in B are reordered 
and switched given a new matrix B sort such that AB' sort ~ I. 
Finally, the A uncertainties are estimated by calculating the 
sample standard deviation of the generated set {B solt }. 

4. STATISTICAL RESULTS 

4.1. Statistical Description of Burst Parameters 

In Table H we present basic statistics of the burst parame- 
ters based on our sample of 197 GRBs. The first two columns 
show the median and the quartile deviation (QD). Since most 
parameters have approximately log-normal distributions, they 
are robust estimators of the typical values and dispersions. 
The ranges are given in the next two columns showing the 
minimum and maximum values. For later reference, we show 
in the fifth column the standard deviations a of the logarithm 
of each given parameter. We tested the observed distributions 
for log-normality using the Kolmogorov-Smirnov (K-S) test 
which is a st andard statistical test to compare any two distri- 
butions (e.g. JPress et all!992l) . and the Geary (G) test which 
is specifically design to test the normality distribution case 
dDevorell982l) . The last two columns list the estimated prob- 
abilities. Note that since we derived from the sample the mean 
and dispersion parameters of the analytical log-normal distri- 
bution, the K-S test probabilities had to be estimated from 
Monte Carlo simulations. Variables like ?9o, r, and Sp seem 
well described by the log-normal distribution, while for Ti ag 
and a it is likely a poor approximation. Nevertheless, in all 
cases it provides better fits than the normal distribution. 
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Some of these GRBs parameters have been extensively 
studied in the literature (see references in § 12.21 , This is cer- 
tainly the case for 790, 7go, T, E p y, and a. The ranges and dis- 
tributions observed in our GRB sample are consistent with the 
ones derived in these earlier works for long GRBs. The only 
exception is the bias towards larger fiuences T due to our se- 
lection criteria. For the new (Sp and TZE p ^) or the less studied 
(t, V, and ri ag ) burst parameters, we show in Figure^their re- 
spective probability density functions (PDFs) and cumulative 
distribution functions (CDFs). The latter are compared with 
the best-fit log-normal distribution (dashed lines). 

4.2. Correlation Matrix 

In Table we show the correlation matrix for our set of 10 
variables. For a sample of 197 bursts, the correlation coeffi- 
cients R = 0.18,0.23,0.27 have a probability of chance occur- 
rence p = 0.01,0.001,0.0001 respectively. We estimated un- 
certainties using the bootstrap method. The correlation matrix 
shows few strong correlations, and many weak but significant. 
As always, it should be remembered that correlation coeffi- 
cients are often sensitive to observational biases and selection 
effects. When the actual correlation involves more than two 
variables forming, e.g., a planar distribution of the data points 
in a multidimensional space, the result of a bivariate analysis 
depends on how these data points are distributed in that plane. 
In contrast, the normal to that plane is rather insensitive to 
the actual distribution of points so long as the thickness of the 
plane is much smaller than its extent. In spite of these lim- 
itations we include here a few comments on the correlation 
coefficients as presented in Table |2] This is done in order to 
compare our results with those of previous work. 

The strongest correlation coefficient R is between 7go and r 
with /?[iogT 50 ,iogr] = 0.87. Naturally, the correlations between 
these two and the rest of the parameters are very similar; how- 
ever, it is important to notice the differences that appear in the 
relation with the other two temporal parameters 790 and V, 
with 750 showing significantly stronger correlations in both 
cases. The fluence T strongly correlates with the temporal pa- 
rameters. These fluence correlations have been reported ear- 
lier in studies of the BATSE catalog, and they are only slightly 
increased here by our sample selection. With our selection cri- 
teria, the lower threshold of the fluence depends on the dura- 
tion and, hence, there is a detection bias against weak and long 
bursts. It is difficult at this stage to determine how much of 
the correlation is produced by this bias; however, in an anal- 
ysis of the fluence truncation effects Lee & Petrosian|jl997fc 
concluded that it is mainl y intrinsic (see also [Mukheri ee et al.l 
1998: lBalazs et al.ll2003t) . The observed correlation between 
T and E p ^ is equally well-k nown, and it was also shown to be 
primarily of intrinsic origin ( Mal lozzi et al.ll995tlLlovd et alJ 
2000). A least-squares fit gives E v ^ oc ^■ 29±Q 03 in ful l agree- 
ment with the index 0.28 ±0.04 found bv lLlovd etafl (2000). 

The newly defined parameters, the emission symmetry Sjr 
and the peak energy ratio TZE p ^, show a significant cor- 
relation. Such a correlation is expected in simple bursts 
with FRED-like (fast-rise slow-decay) structure due to the 
observed hard to soft e volution during the decay phase 
(Borgonov o & Rvdell2001l) . Hence, this correlation indicates 
that the more FRED-like is the light curve the better is the 
hardness-intensity correlation. These new parameters also 
show an anti- corre lation with the variability V, that is further 
discussed in § 15.11 There seems to be only a weak correlation 
between the two spectral-evolution parameters TZE p ^ and ri ag . 
The riag also correlates weakly with (-)£pk and a, and it ap- 



pears to be independent of the other parameters. iNorris et 
(2000) found, for a small sample of GRBs with known red- 
shift, an anti-correlation between peak luminosity L and lag 
r/ ag calculated in the local frame. However, as the authors 
noted, in the observer's frame there is neither an apparent 
correlation between peak flux F and lag Ti ag (we calculated 
for our sample R{\ gF,\o%n„] = 0.004), nor between fluence and 
lag (see Table |2j. In order to verify that our result is consis- 
tent with that of Norris et al. (2000J), we used the observed 
sample of lags T\. dg together with redshifts for 37 pie-Swift 
GRBs and the bootstrap method to randomly draw surrogate 
data. We then generated realizations of the lag-flux data sam- 
ple by applying th e power-law relation derived for L(t') by 
Norri s etal] (J2000), adding noise to reproduce the observed 
dispersion. Taking into account the sample flux range, no sig- 
nificant correlation was obtained in agreement with the obser- 
vations. We find a weak correlation between T and V, which 
is probably a consequence of the c orrelation between T and 
750 (as we will show later in § 153 and most likely has no re- 
lation with the luminosity-variability correlation claimed by 
iReichart etal] (2001) since for our sample /?[i og F,iogV] = 0.01. 
We also f ind an anti-cor relation between V and ri ag , as re- 
ported by Scha efer et alJ d2001h . The low-energy power-law 
index a shows only weak correlations. The known anti- 
correlation between a and E p ^, although here hardly signifi- 
cant considering the uncertain ties, has bee n explained as a fit- 
ting a rtifact (see, e.g., iPreece et"ai]ll998tlLlovd & Petrosianl 
2000). 

4.3. Principal Components 

From the correlation matrix shown in Table [2] we derive 
the matrix of coefficients A (defined in eq. Q) that describes 
the linear transformation of the initial variables into the ba- 
sis of PCs. In Table |3 the z'-th data row corresponds to the 
PCj indicated in the first column. The coefficients a,j from 
equation (|2j that determine the PCs are given in columns 2- 
11, where each column is associated with one of the variables 
in Table [2] The last two columns show the percentage of to- 
tal variation, Var, and the cumulative percentage of variation, 
CVar. 

To estimate the corresponding uncertainties we used the im- 
plementation of the bootstrap method described in §|3] which 
takes care of possible re-ordering and sign switching of the 
PCs. Nevertheless, some degree of "mixing" cannot be pre- 
vented when two eigenvalues are equal to within uncertain- 
ties. That is the case between PCs and PC(„ where the errors 
are consequently overestimated (i.e., the shown uncertainties 
should be taken as upper limits). On the other hand, PC] 
and PC10 are examples of "stable", well-determined principal 
components. Additional degeneracy may exist in subspaces 
associated with subsets of variables that are closely interre- 
lated (e.g., the subset of temporal variables). Unusually large 
uncertainties in some coefficients within a PC are likely due 
to this problem. 

As we explained in § [3] one of the PCA applications is the 
derivation of a reduced set of variables to de scribe the data . 
Adopting the conservative criterion given by Jolliffe (2002), 
for a set of n variables, the number m of PCs needed to re- 
produce the observe variation of the data is given by those 
PCs with more than (70/«)% percentage of variation. In the 
present case n = 10, so looking at the percentage of varia- 
tion Var in Tabled a rejection level of 7% gives m = 5. This 
suggests retaining the first 5 PCs, although the last two are 
marginal. The implications of these PCs are discussed in more 
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FIG. 1 . — Probability density functions (PDFs) and cumulative distribution functions (CDFs) for five of the burst parameters included in our analysis (i.e., the 
autocorrelation function (ACF) half-width at half-maximum r, the variability V, the emission symmetry Sjr, the time lag tw, and the peak energy ratio TZE^). 
In the CDFs we show for comparison both the observed distribution (solid lines) and the analytical log-normal distribution that best describes it (dashed lines). 
Very good fits are obtained for the logarithms of r and Syr (see TablefTlfor the results of comparative statistical tests). Note that the newly define parameters Sjr 
and TZEpk are in most cases (but not always) larger than unity implying, according with their definitions, that for most bursts there is an overall flux decay and 
hard-to-soft evolution. 



detail later in § 15.31 Instead we will start with the analysis of 
the PCs with the smallest eigenvalues, which have a more di- 
rect physical interpretation. 



5. ANALYSIS 



All burst parameters are subjected to some degree of ob- 
servational bias. This can have many causes like the sensi- 
tivity limit of the instrument, its energy band, the temporal 
resolution of the light curves, etc. As mentioned above such 
biases together with the selection criteria can sometimes sig- 
nificantly distort the observed relation between two param- 
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FIG. 2. — Relation between four temporal variables: duration time Tgo, emission time 750, half-width at half-maximum of the autocorrelation function r, and 
high frequency variability V (bluish and reddish dots indicate high and low variability values, respectively). The data points are distributed in an angular sector 
where the lower limit is given approximately by the identity relation (dashed line). A power-law fit of the data (solid line) gives an index 0.60 ± 0.03 in agreement 
with equation |4j. 



eters and generate false correlations. In contrast, the direc- 
tions of the PCs with the smallest eigenvalues (see § |3j are 
more robust against systematic errors. Even if a threshold 
in one of the parameters cuts off part of the sub-space or- 
thogonal to such a PC, or a bias affects the way points are 
distributed within that sub-space, the orientation of the sub- 
space is unlikely to be substantially affected. Hence, analyti- 
cal constraints imposed on the variations of parameters from 
PCs with small eigenvalues are expected to be rather insensi- 
tive to biases and selection effects. 

5.1. Tgo-Tso-T Correlation 

The last PC has associated with it a very small percentage 
of the total variation, with Var ~ 0.5%. For this reason PC\o 
represents a nearly constant relation between the few variables 
with significantly large loadings, i.e., the temporal variables 
790, 750, and t. Taking into account that the coefficients in 
Table |3] give the linear combination of the logarithms of the 
original variables, normalized by their standard deviations a, 
it follows 

0.36 0.79 045 

T g ° ( f T 50 0Ai r^9 ~ const. , (3) 

where the er-values were taken from Table ^ Then equa- 
tion Q can be approximately rewritten as 

(T 5 o/r 9 o)cx(r/r 9 o) - 6 , (4) 
so that the duration of the burst can be used as normalization. 
The main thing to notice from equation (0} is the indication 
of a self-similar structure of the GRB light curves. This is a 
non-trivial property; two of the exponents in equation (0 are 
independent and, hence, it is in general not possible to arrange 
three PCA parameters in this way. It is also clear that this is an 
intrinsic property, since the ratios in equation (0} do not de- 
pend on redshift. We should also remark that the measure r in 



long burst light curves is at least several times smaller than the 
time duration Tgo- Therefore, the r time-scale is in most cases 
well sampled, and it is not significantly affected by the limited 
burst duration. Furthermore, it directly relates to the observed 
low-frequency break in the typical power-law behavior fou nd 
in power density spectra analysis ( Belo borodov et aiT2 000). 

In order to investigate the generality of the self-similarity in 
equation we have subdivided our sample, for examples, in 
halves of high and low values of 790. The same self-similar 
structure is found in the subsamples, and the same relation 
between the temporal variables is derived when performing 
separated PCAs, although with lower statistical significance 
due to the smaller number of sources in each subset. 

In Figure 12 we show a scatter plot of equation (0}. Con- 
sidering that the coefficient associated with the variability V 
in fCio has a small but non-negligible value (see Table [5), 
we indicate also in the figure the degree of variability V us- 
ing a color scale, relating in this way all four temporal vari- 
ables. The data points lie within a roughly angular sector 
where the lower limit is given approximately by the identity 
relation 750 ~ r (dashed line). When fitting the data points 
using a power-law model (shown in solid line) we obtain an 
index 0.60 ± 0.03 in agreement with equation (|4j. It is seen 
that bursts with low values of V (reddish points) have a distri- 
bution which differs significantly from those with high values 
(bluish points). Bursts with large values of V are mainly found 
well above the identity line and close to the upper-right ver- 
tex, while those with low values of V tend to lie closer to the 
identity line and at some distance from the vertex. Since the 
value of V is a measure of the complexity of the short time- 
scale structure in a given burst, this indicates that the distribu- 
tion of V-values in Figure|2]could be d ue to a vary ing number 
of pulses (see also the discussion by iLi & P aczvnskii l2006l 
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FIG. 3. — Synthetic data points based on artificial burst light curves. A 
simple toy model is used to generate light curves with a random number of 
exponential pulses (see text for details). The range and distribution of the data 
points are well matched with those in Figurel2lfor « max = 25 and a power-law 
fit (solid line) gives an index 0.58 ± 0.03 in agreement with equation |4j 
within uncertainties. The color scale indicates the number of pulses. Bursts 
with few pulses (red points) lie close to the identity relation (dashed line), 
while bursts with a larger number of pulses lie further from this line (orange- 
green points). As overlap becomes increasingly important (green-blue points) 
in bursts with a large number of pulses, they move towards the upper-right 
vertex. In Figure|5]the regions of high and low variability V follow the same 
behavior pattern, indicating that V is an indirect measure of the number of 
pulses. 



on the connection between the value of V and the number of 
pulses). 

In order to test this idea we proceeded to reproduce the ob- 
servations using synthetic burst light curves based on a toy 
model that generates a random number of exponential pulses. 
We found that the qualitative behavior does not depend on 
the pulse shape. To minimize the number of free parameters 
we fixed the pulse amplitudes and the burst duration Tgo, the 
latter motivated by the observed self-similar structure of the 
bursts. We adopted a normal distribution for the decay char- 
acteristic times fd- The parameters of the fd distribution (i.e., 
its mean value and dispersion) are basically set by those of the 
observed r-distribution since fd ~ r. Furthermore, a uniform 
probability distribution (7[l,n max ] was assumed for the num- 
ber of pulses. Values for fd and the number of pulses were 
chosen randomly from their assumed distributions. Bursts 
were then generated by placing the pulses randomly within 
the time span Tgo. We adjusted the only free parameter n max in 
order to approximately reproduce the range of values, means, 
and dispersions observed in the 750 ordinate, and to obtain a 
power-law fit with index ~ 0.6 as in equation 

The best match was found using a pulse distribution with 
"max ~ 25. With this simple model we were able to repro- 
duce all the main features observed in Figure |2] An exam- 
ple of a good match between the outcome of our simulations 
and the observations is shown in Figure [3] where the color 
scale now corresponds to the number of pulses. Single-pulse 
bursts are found close to the identity line because in those 
cases 750 ~ t ~ fd- Bursts whose light curves have a few 
well-separated pulses lie above the line since for uncorrelated 
random pulses r will be approximately the same as for a sin- 
gle pulse, while the emission time 750 will increase (red and 
orange points). The chance of pulse overlap grows with the 
number of pulses, and as they start to overlap a shift towards 
larger r-values occurs and, hence, both r and 750 increase. 
Finally, bursts showing many heavily overlapping pulses re- 
semble a single broad pulse and these will populate the upper- 
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FIG. 4. — Analysis of the data-point distribution in Figure |2]taken the ra- 
dial distance R from the upper-right vertex as reference. On the log-linear 
histogram representation the density decrease is approximately linear (solid 
lines). Bursts with heavily overlapping pulses lie near the vertex at small 
radii, where the relation between an d T "saturates" (see text). The dis- 
tribution for the synthetic data in Figure|3]is shown for comparison (dashed 
lines). Based on the K-S test the differences between the two distributions 
are not significant. 



right vertex (green and blue points). In Figure|2]the regions 
of high and low variability, V, follow the same behavior pat- 
tern, indicating that V is an indirect measure of the number of 
pulses. The bursts with low values of V at the upper-right ver- 
tex would then correspond to light curves with a large number 
of pulses so that the pulse overlap is large enough to resem- 
ble a single pulse. A complementary way of comparing ob- 
servations with simulations is used in Figure [4] which shows 
the radial distribution of points taken the upper-right vertex as 
reference. 

The main features of this toy model rely on the self-similar 
characteristics of the statistical properties of r and T50 (i.e., 
their mean values and dispersions scale with 790). This im- 
plies that significant overlap of pulses start to occur when 
the number of pulses is n , where n oc (r/Tgn) -1 / 2 (assum- 
ing « 3> 1). The value of 750 at this point is, roughly, n r, 
or T 50 (xt 1 / 2 . This power-law index in the relation between 
750 and r is smaller than the observed one (ss 0.6), indicating 
that overlap of pulses is more significant for larger values of 
r/TgQ. In the toy model this is achieved by assuming « max to 
be independent of r/Tgo- This is a non-trivial aspect of the 
toy model, since this independency together with the actual 
value of n max are both needed in order to account for the two 
independent characteristics of the observed relation between 
Tgo and r, namely, the power-law index and the radial distri- 
bution shown in Figure 0] 

5.2. Epk-J-'T Correlation 

The ninth PC has also a small associated variation (Var ~ 
2.1%), and disregarding small loadings as before, it can be 
interpreted as a fairly constant relation between T, £ p k, and 
r. Since the correlations of 7so with the first two variables 
are the same as the corresponding ones with r, some degener- 
acy between the two temporal variables could have been ex- 
pected. However, the associated uncertainties are not partic- 
ularly large which indicates that, in the context of this 10- 
dimensional multivariate analysis, it is r that best maximizes 
the variance of PCg. On the other hand, doing the PC A with 
only three variables gives almost identical results using either 
r or T50. The relation between the three original variables can 
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be approximately written as 

£ pk oc ^°- 8 T - () - 6 , (5) 

where the reported exponents have been derived from a PCA 
restricted to the three-variable subset. They are a more ade- 
quate reference for the numerical simulation that will be pre- 
sented below, which is limited to the same subspace. In any 
case, the index values in equation (0 are equal within uncer- 
tainties to the ones that are obtained from Tables|2]and[3] 

The fact that the strength of the correlation between E p ^ 
and T is increased by an additional temporal variable is par- 
ticularly interesting, since an analogous behavior has been 
obser ved in the study of the corresponding local frame vari- 
ables. lAmati et all (|2002) found that in the local frame E' k oc 

£'°q 2±0 06 , where E lso is the i sotropic equivalent gam ma-ray 
energy. It was later shown by Ghirland a et alJ fc00 4) that the 
scatter around this relation derived mainly from different val- 
ues of the break time, as measured in the afterglow light 
curves. This suggests a relation between characteristic time- 
scales in the prompt emission (r) and in the afterglow (f^). An 
additional thing to note from equation (0 is the exponent of 
T, which differs markedly from the one deduce above from 
the bivariate correlation, namely E v \ oc jr0.29±o.03 j n ^ s jjj us . 
trates the fact that projection effects can strongly influence the 
deduced relation between observables. 

The observed relation between E p ^, T, and r in equation 
as well as the observed variance are affected both by the un- 
known redshifts and the sample selection criteria. Hence, the 
usefulness of equation depends critically on the possibil- 
ity to account for these effects. To this end, we assume a local 
frame relation of the form 

^ k cx£f so r' 6 , (6) 

for which we like to determine values of a and b consistent 
with observations. Since the inversion problem has no unique 
solution under the present conditions, the use of Monte Carlo 
methods together with a forward-folding technique seems to 
be the best way to constrain the {a,b} parameter space. Re- 
stricting the PCA to a 3-dimensional subset, the third PC gives 
the relation shown in equation 0, but with a higher associ- 
ated eigenvalue of 7.7% since now this percentage is taken 
over a smaller total variation (i.e., from Table y] approxi- 
mately 10/3 times 2.1%). We will assume that all source and 
observer's frame variables have log-normal distributions and 
therefore are fully determined by their means and standard 
deviations. 

In principle, the PCA can be used to simulate any correla- 
tion matrix for the synthetic data, generating first random re- 
alizations of the uncorrelated PCs and then transform them to 
the original variables using the matrix of eigenvectors. How- 
ever, in this way it is difficult to explore systematically the pa- 
rameter space {a,b}, since these parameters can not be used 
as starting conditions. Initially, we tried this general approach 
but only to verify that E pk and r' can be assumed uncorrelated, 
and a simpler procedure was then used for our problem. 

We proceeded instead by estimating values for the param- 
eters of the probability distributions of E pk and r' and used 
these to generate random values. The E^o-value is then de- 
rived using equation for a given (a,-, pair. In this way we 
are exploiting also the fact that Ej S0 has the largest dispersion. 
It is assumed that the rest-frame relation in equation has 
no intrinsic variance, i.e., that £i so , E' pk , and r' lie strictly in 
a plane. Since the observed variance in equation is a con- 



volution of the intrinsic variance and that introduced by un- 
known redshift and selection criteria, the derived variance is a 
lower limit. Hence, acceptable values for {a,b} are those that 
in addition to reproducing equation 0, also have a variance 
not larger than the observed one. Next, we convert the data to 
the observer's frame. To do so, we estimate the redshift dis- 
tribution from the pie-Swift sample of known redshifts. The 
conversion T = Ej S0 ( 1 +z)/ (47rd|) involves the luminosity dis- 
tance diXz), which depends on cosmological parameters. We 
adopt current standard values with a Hubble constant Ho = 65 
km s" 1 Mpc" 1 , a matter density tt m = 0.3, and a vacuum en- 
ergy density 51a = 0.7. We also consider E[ so , E' k , and r' 
redshift independent (i.e., evolution effects are neglected). 

The simulated data point (JF ,E$.,t) ^ thus obtained is then 
accepted as an observation or rejected according to the obser- 
vational constraint s of our sample data and the sample selec- 
tion criteria (see § 12.11 . The BATSE energy window is much 
wider than the observed E v y distribution, so the selection is 
not very sensitive to its effective limits. However, the treat- 
ment of the flux selection criterium is less straightforward, 
since a value for the flux is not directly available in our simu- 
lations. Hence, the flux selection criterium has to be translated 
into a criterium relevant for the parameters in the simulations. 
In a bivariate correlation analysis between the flux F\ s and 
each one of our variables, only the fluence T seems to have 
a significant correlation with /?pogfi s ,iog^] = 0.56. However, 
a multivariate PCA involving the flux and the three parame- 
ters under study reveals a well-defined thin correlation-plane 
in the {Fi s ,F, r} parameter subspace (with associated eigen- 
value 2%). Note that we do observe a correlation between 
T and r (Table 0, but we found a negligibly small lvalue 
between F\ s and r. This illustrates that in general, when deal- 
ing with real experimental data, the measure of correlation R 
does not necessarily show the transitive property. It is also 
a good example of the robustness of the multivariate method 
discussed before in §0 The intersection of the {Fi s ,JF ; r} 
correlation-plane with the boundary-plane given by the selec- 
tion criterion F\ s > 4 photons cm 2 s _1 determines the threshold 
of acceptance. This can be approximately modelled using the 
inequality log J 7 > (-5.5 + logr). 

After the acceptance test, the process is repeated until the 
number of data points equals the size of the sample (197). 
Then, means and standard deviations are calculated and com- 
pared with the observed values and the initial parameters of 
the local distributions of variables are corrected accordingly in 
an iterative cycle until the calculated values agree with the ob- 
served ones. By this method our estimation of the E{ &0 disper- 
sion is consistent with th e ones based on the few GRBs with 
known redshifts (iFrail et al.l200lHGhirlanda et all2004l) . 

Once the appropriate distribution parameters for the local 
variables are estimated, it is straightforward to generate many 
fiducial (FjEpbjT) sets. We then apply the PCA to each set 
and calculate means and standard deviations of the eigenvec- 
tors and eigenvalues. In particular, the values associated with 
the third PC (i.e., the last PC in this three variables PCA) 
are compared to the observed relation (eq. 0). If they agree 
within uncertainties, then the initial values (a,-, bi) are accepted 
as compatible with the observations. Values of (a,-,£>,) consis- 
tent with equation are shown in Figure 0together with the 
associated eigenvalues (variances). 

As mentioned above, valid eigenvalues should not be larger 
than the observed value (7.7 ±0.7)% within uncertainties (ap- 
proximately 10% of the reported values). As can be seen from 
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FIG. 5. — Parameter space of the power-law indices {a,b} in the relation 
£p k oc E" so T ,b . The dots indicate values compatible with the observations 
(i.e., eq.[|5)), and they are labelled with the corresponding PCA eigenval- 
ues that represent the percentage of the total variation. The best match with 
the observed eigenvalue (7.7 ±0.7)% is (8.1 ±0.8)% which corresponds to 
(a,b) = (0.5,-0.3). 



Figure[5] in all cases the calculated eigenvalues are equal to or 
larger than the observed ones, which implies: 1) Acceptable 
values of {a,b} lie within a small region of parameter space, 
and 2) the intrinsic variance in the relation between E[ so , E^ k , 
and t' is considerably smaller than the observed one, which is 
then mainly due to unknown redshifts and the sample selec- 
tion criterium. 

The above conclusions regarding the intrinsic properties of 
the relation in equation (jfji are quite restrictive. However, they 
do rely on the validity of our treatment of unknown redshifts 
and sample selection criteria. We now discuss these aspects 
in turn. It could be argued that the set of pie-Swift GRBs 
with known redshifts z used above is not representative of the 
whole BATSE sample of long GRBs, since the determination 
of z always involved the more precise localization provided 
by other instruments with different trigger criteria and lower 
sensitivities. However, it has been shown by statistical anal- 
ysis of the distributions of physical parameters that there are 
no significant differe nces between that s et and the set of long 
bright BATSE bursts ( Bos niak et alJ2005h . which are the ones 
considered in this work. In order to directly investigate the in- 
fluence of the redshift distribution on our results, a broader 
redshift distribution that peaks at z ~ 2 (i.e., about twice the z 
mean value and dispersion observed by pre-Swift instruments) 
was assumed. This gave very similar results to those pre- 
sented in Figure[5] although with somewhat larger uncertain- 
ties associated with the eigenvalues. Hence, the main con- 
clusions of our analysis are not very sensitive to the assumed 
redshift distribution. Furthermore, among the different pos- 
sible sample biases, we found the empirical modeling of the 
selection effect created by the correlation plane {F\ s , J 7 , t} to 
be the only relevant one. However, for reasonable choices of 
the threshold J-(j), we found small variations in the index b 
and, in particular, a value significantly larger than zero. We 
therefore conclude that b = can be rejected and, hence, that 
a significant fraction of the scatter in the well-known Ei so -E' k 
relation must be due to variations in the value of r'. Our best- 
estimate values of the power-law indices in equation (|6) are 
a = 0.50 ± 0.05 and b = -0.30 ±0.05. N ote that the £ iso index 
is in agreement with Ama ti et all (2002). 

5.3. Correlation Structure 



While the last PCs lead to direct functional dependencies 
between various variables, the first few PCs instead reflect 
global properties in that sub-sets of variables are identified 
which are more strongly interlinked with each other than with 
the rest. If all the variables were equally contributing to the 
first PC, it would have approximately equal coefficients or 
loadings, i.e., considering the normalization for n = 10, each 
PC\ coefficient would be « 0.32. Instead, for PC\ the load- 
ings are dominated by the temporal parameters ?9o, %o, t, and 
to a lesser extent V. In PC2 we see roughly the complemen- 
tary pattern, with significant loadings for the spectral group of 
variables, in particular E p ^ and to a lesser extent IZEpk and a, 
but also Sjr. The only parameter with significant loading in 
both of these PCs is T. The vector PC3 is mainly driven by 
the two spectral evolution variables T\- dg and TZE p ^. The vari- 
able V seems to contribute non-negligible loadings to all three 
PCs, although with somewhat lower significance for PC2 and 
PC3 due to the larger error bars. One may also note that the 
loadings for V occur with opposite signs to those of and 
7\L£pk, indicating an anti-correlation. Since Sp- is relate d to 
the FRED-like appearance of the burst, it was argued in § 14.21 
that the degree of correlation between and TZE p ^ is a mea- 
sure of the prominence of the hardness-int ensit y correlation 
in the burst. In the toy-model discussed in § 15.11 V is directly 
related to the number of pulses in the burst. Hence, if the 
hardness-intensity correlation as well as the FRED-like ap- 
pearance are associated with the individual pulses rather than 
the burst itself, the fact that V anti-correlates with both S?- 
and IZEpk could be attributed to the expected weakening of 
both these characteristics as the number of pulses increases. 

One of the objectives of empirical methods like the mul- 
tivariate correlation analysis is to reveal the multiple depen- 
dencies between the observed parameters. The PCs with 
the smallest and largest eigenvalues give complementary in- 
formation regarding the statistical properties of the sample. 
However, such an analysis has many limitations; for example, 
one cannot tell whether a correlation that is found indicates a 
direct causal connection or just an association through some 
variable, not necessarily included in the analysis. In order to 
give a somewhat different view of the statistical properties of 
our sample of GRBs we tried various other methods. The sim- 
plest approach was to directly inspect the correlation matrix 
(Table |2j, ranking the strongest correlations for each param- 
eter. This method is illustrated in Figure |6ji. If a given cor- 
relation is ranked high for two parameters it would indicate 
a close connection between them. However, this reciprocity 
is rare when dealing with a large number of variables and in- 
stead it is more likely to find, like in our case, subsets that are 
interlinked. The problem then is that one cannot distinguish 
between indirect and weak dependencies. 

An intermediary method, which combines the ranking of 
correlations and the PCA, is to use the PCA to find the cor- 
relation between a given parameter and any number of the 
remaining parameters. The simplest version is to find three- 
parameter correlations and to derive relations like equation (|5} 
where one of the parameters appears as a product of power- 
law functions of the other two. Hence, two parameters are 
used to form an estimator of the third and the best model is 
defined as that giving the largest correlation coefficient be- 
tween predicted and actual value. On several numerical ex- 
periments with multiple intercorrelated random variables, this 
scheme was found useful to reveal which variables were more 
directly connected. A variety of functional dependencies were 
investigated, from simple "chain" structures (where one vari- 
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FIG. 6. — Schematic representation of the correlation structure for our set of 10 GRB parameters. Reciprocal relations suggest a close proximity between 
parameters (solid lines), while non-reciprocal relations indicate a correlation of lower significance (dotted lines with single arrows), (a) Each parameter is 
pointed at by the parameter(s) showing highest correlation with it, as given by Table l2l (b) We apply the PCA to all possible triad subsets and find for each 
parameter which other two are the best to construct an estimator. In the context of this analysis, for the fluence T equally good estimators can be derived using 
£pk and either r or 750- The numbers indicate the correlation coefficient R of each parameter with its corresponding two-parameters estimator. 



able drives another that in turns drives another, and so forth) 
to more complex structures generated by multiple crossed as- 
signments between the variables. Random noise was added 
at every instance to weaken the correlations. In all cases the 
PCA based method consistently produced better results than 
the direct reading of the correlation matrix (i.e., the two di- 
mensional case). Under high noise levels the test produced 
less but rarely false identifications. 

For each of the GRB parameters in our set of 10, we con- 
sider all possible three dimensional subsets and then calculate 
their PCs. In each case we derive a parameter estimator from 
the last PC (i.e., the one with the smallest eigenvalue) and cal- 
culate the corresponding correlation with the actual parame- 
ter. In this way we find for every parameter which other two 
produce the best estimator. The analysis is summarized by 
the graph shown in Figure[6j?. Reciprocal relations are taken 
as evidence for a close association and the parameters are ar- 
ranged accordingly. Non-reciprocal relations (also shown) in- 
dicate associations at a significantly lower level of confidence. 
The first to notice from the two graphs in Figure [6] is that the 
use of the PCA-method reveals considerably more structure 
than does the ranking of correlations. The two distinct subsets 
of parameters that are apparent in Figure[6j? directly relate to 
the first three PCs. However, Figure [6{? provides information 
in addition to that contained in the first PCs, since it explic- 
itly shows the relationship between the variables within each 
of the two subsets and also, to some extent, the connection 
between the subsets. The subset of temporal parameters is 
equivalent to PC\ , while the parameters in PC2 and PC3 seem 
to be joined in the other subset indicating a close connection 
between them. It is seen from Table [3] that the eigenvalues 
for PC2 and PC3 are rather similar and have values roughly 
half of that for PC\ . This suggests that although the spectral 
evolution is related to the overall spectral properties, it has 
its own distinctive features. Two additional things to noticed 
from Figure [6j> is that Sjr belongs to the second subset of 
parameters and that V appears to have no strong relation to 
either of the two distinct subsets. This argues that not only 
the FRED-like structure but also the basic temporal and spec- 
tral properties are associated with the individual pulses. The 
overall properties of the burst is then determined mainly by 
the number of pulses. 

6. DISCUSSION 



The statistical properties of a sample of GRBs have been 
analyzed with the use of ten parameters chosen to bring forth 
the main characteristics as observed by BATSE. In a multi- 
variate analysis it was found that ~ 70% of the total vari- 
ation in parameter values is driven by only three principal 
components. Half of that is due to temporal variations while 
the other half is roughly divided equally between two prin- 
cipal components describing, respectively, the average spec- 
tral properties and spectral evolution. A more refined analysis 
reveals that the parameters in the latter two principal com- 
ponents are closely related. Together with the similarity of 
the variance (eigenvalue) of these principal components, this 
makes it likely that the two principal components needed to 
describe the spectral variations are not uniquely defined by 
the present analysis. Hence, the conclusion is that there ex- 
ist two distinct sub-groups of inter-related parameters, which 
describe the temporal and spectral properties, respectively, of 
GRBs. One may note from Figure[(J> that the only parameter 
showing substantial connection to both sub-groups appears to 
be the fluence, i.e., its value is determined jointly by the two 
sub-groups. 

Complementary information is provided by the principal 
components with the smallest variance. Since these direc- 
tions in parameter space are defined by their narrow distribu- 
tions of parameter values, they directly provide analytical re- 
lations between the parameters involved. It was found that the 
parameters in the temporal sub-group discussed above also 
obey such a well defined relation. It is noteworthy that this 
relation is self-similar in the sense that the mean values and 
dispersions of both the emission time (Tso) and half-width of 
the autocorrelation function (r) scale with the duration of the 
burst (7go). Since this relation is unaffected by the unknown 
redshifts, it is a characteristic intrinsic to the light curves of 
GRBs. It is seen from Figure[6j> that the parameter V, which 
measures the small time-scale variability, is the only one (pos- 
sibly also a) without any obvious c onne ction to either of the 
two sub-groups. It was shown in § !5.1l that the observed re- 
lation between 750 and r finds a ready explanation in a toy 
model where the value of V is determined by the number of 
pulses in a burst. This implies that the scaling of 750 with Tgo 
would be the same as that for t and, hence, the fundamental 
scaling is that between r and 790. 

Burst light curves have very di verse morphologies and 
it has been suggested (see, e.g., iGolenetskii et alJ Il983t 
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iBoreonovo & Rvdel2 001 ) that they are a composite of pulses, 
with each pulse representing a single physical event (e.g., in 
the internal shock scenario, collisions between shells). In the 
toy model, the number of pulses in a burst is measured by the 
value of V. Independent support for such an interpretation 
comes from the prominence of the hardness-intensity corre- 
lation. It is well-known that this correlation is most apparent 
in FRED-like bursts. This association is most clearly seen 
in our data when these are represented as in Figure [6j? (i.e., 
the close relation between Syr and TZE V ^). The effects of V 
can best be seen in the PCs with the largest variances; in all 
three, V anti-correlates with both Sj? and 1ZE V ^, indicating 
that large values of V tend to smear out both the FRED-like 
structure and the hardness-intensity correlation. Hence, a sim- 
ple interpretation of the data is that not only the self-similarity 
of the light curves but also the hardness-intensity correlation 
and the FRED-like structure are properties associated with the 
individual pulses and that the overall properties of a burst is 
determined by the number of pulses. 

The correlation between fluence (J 7 ) and peak energy (£pk) 
was one of the first to be well established for GRBs. Subse- 
quent analysis indicated that the scatter around this relation 
could be reduced if a time variable were included, although 
it remained unclear which time variable gave the best result. 
Part of the problem, when trying to make more quantitative 
estimates, was the lack of redshifts. Wit h the use of a small 
sample of GRBs with known redshifts, lAmati et"ail ll2002l) 
found the rest-frame equival ent of the ■7 r -fi' n i c -relat ion (i.e., 
the Ziiso-iipk-relation), while Ghirlanda et al. (2004) showed 
that the scatterer around this relation could be substantially 
reduced if account were taken of the rest-frame break-time 
(tj) in the light curv e of the ensuing afterglow (see also 
Liang & Zhandl2005h . It is interesting to note that the func- 
tional dependence between E p ^ and E mo is the same in the 
Amati and Ghirlanda relations. This is due to the fact that the 
values of t' h and E' k are uncorrected. The reduction of the 
scatter in the Amati-relation with the use of a time-scale as- 
sociated with the afterglow rather than the prompt emission 
indicates a close connection between these two phases. 

Although the lack of redshifts prevents a determination of 
the relation between the values of E lso and E pk for individ- 
ual GRBs observed by BATSE, we have developed a method 
by which this relation can be investigated statistically for the 
sample as a whole from the observed E^-T-t correlation. 
By matching both the mean values and the dispersions of 
these observables, we show that the Amati-relation can be 
derived from our sample of GRBs. Furthermore, the scatter 
around this relation correlates with the value of t', the rest- 
frame value of the half-width of the a utocorrelation function . 
This is analogous to the findings of Ghirla nda et al.l (2004) 
that the scatter in the Amati-relation correlates with t^. The 
observed scatter in the E^-T-t correlation is consistent with 
that expected from the unknown redshifts, indicating a small 
intrinsic scatter. During the prompt phase of a GRB, it is 
likely that the observer sees only a small portion of the jet; 
hence, r' measures a local property and its value could, in 
principle, vary over the jet (i.e., for different lines of sight). 
The value of f£, on the other hand, is due to the overall prop- 
erties of the jet. The fact that both r' and reduce the scatter 
in the Amati-relation indicates that at least some of the local 
p roperties do not vary m uch over the jet. 

INakar & Pirar] ( 120051) and lBand & Preecd ( 120051) have ar- 
gued that the Amati-relation is not consistent with the BATSE 



data and may be the result of a selection effect. However, as 
shown here, the Amati-relation can, in fact, be derived from 
our sample of GRBs selected from the BATSE data. This is a 
complete, flux limited sample, which suggests that the Amati- 
relation is valid for the intrinsically most luminous GRBs. 
These apparently contradicting results may be caused by a 
double bias against observing afterglows with large values of 
tj: ( 1 ) Due to the anti -correlation between t' h and E lm found by 
Ghirlanda et al. (2004), afterglows with large values of are 
intrinsically weak. (2) The decline of the light curve with time 
makes it increasingly hard to establish a value for in those 
afterglows where the break occurs late . He nce, the incon- 
sistenc y found bv INakar & Piranl (|2005) and iBand & Preecel 

( 2005 ) could be due to a larger number of GRBs in the BATSE 
sample than assumed by them for which it is observationally 
hard to measure a break in the light curve of the afterglow. 
The Amati-relation would then constitute an upper envelope 
of a continuous distri bution, a possibility also discussed by 
Nakar & Piran (2005). As a result, associating the break in 
the light curve with the opening angle of a jet, implies an aver- 
age opening angle considerably larger than typically deduced 
fr om observations. 

Firm ani et alJ ll2006l) have analyzed a sample of 15 GRBs 
with known redshifts and found a good correlation between 
peak luminosity (L; so ), E' k and 7^' () (actually, they used 7^' 5 in- 
stead of Zj' ). They argue that this is a correlation independent 
from the one between E^-E^-t^ on the ground that the scat- 
ter in the former correlation increases if is used instead of 
7^o as the third observable. This raises the question whether 
7^' and are both related to a third time-scale, not included 
in the analysis, which could then be the appropriate one to 
use in both of these correlations. If this were the case, the 
apparent independency would have no physical ground but 
due to a too limited analysis. Since we have used peak flux 
to select our sample, this observable (and, hence, L ls0 ) was 
not included in the analysis. However, a few of our results 
have bearing on the relevance of the various time-scales. Al- 
though there is a strong correlation between t and T^q, the full 
10-parameter analysis clearly indicates that r is the preferred 
third observable rather than 750 in the £pk-.F-relation. This is 
in contrast to a limited 3-parameter analysis which makes no 
difference between r and T^. Hence, averaging over (or ex- 
cluding) parameters which contribute only a small amount to a 
given correlation, can smear out the distinction between direct 
and indirect dependencies. Furthermore, before applying our 
method to account for the effects of the unknown redshifts to 
the E v y-T-T correlation, we took some care to establish that 
the observed correlation between E^ and r is consistent with 
being due to the unknown redshifts, i.e., that the re is no intrin- 
sic correlation between E' k and r'. As noted by Firmani et al. 

(2006) , their low value for the sum of residuals Xt could be 
due to a correlation between the errors. In fact, in their sam- 
ple the value of E' k correlates with that for 7^' , indicating that 
an observable more fundamental than 7g' may exist (note that 
the values of E pk and t' h are uncorrected). Given our analy- 
sis of the E^-T-t relation, we suggest that this observable 
co uld be t', which was not among the parameters considered 
by Firmani et al. (2006). It is then possible that 7^' and 
are both related to r' and, hence, that E^o-E^-t^ and Li so - 
E pk -T^ are not independent correlations. In this context, one 
may note that in the toy model used above to account for sev- 
eral of the observed properties of GRBs, r' is the fundamental 
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TABLE 2 
Correlation Matrix 



Parameter 


log r 90 


log 750 


logr 


logV 


log <S^ 


log T kg 


log TZE pk 


log.^ 


log £ p k 


a 


log T 90 


l.±0 


0.78 ±0.03 


0.58 ±0.04 


0.18 ±0.07 


-0.09 ± 0.06 


-0.01 ± 0.06 


-0.15 ±0.06 


0.50 ±0.05 


0.24 ±0.06 


-0.26 ± 0.07 


logTio 


0.78 ±0.03 


l.±0 


0.87 ±0.02 


0.51 ±0.05 


-0.25 ± 0.07 


0.09 ±0.07 


-0.21 ±0.06 


0.61 ±0.04 


0.14 ±0.08 


-0.16 ±0.07 


logT 


0.58 ±0.04 


0.87 ±0.02 


l.±0 


0.40 ±0.06 


-0.24 ±0.07 


0.15 ±0.07 


-0.25 ±0.06 


0.61 ±0.04 


0.14 ±0.08 


-0.12 ±0.07 


logV 


0.18 ±0.07 


0.51 ±0.05 


0.40 ±0.06 


l.±0 


-0.32 ±0.07 


-0.18 ± 0.08 


-0.37 ±0.05 


0.33 ±0.06 


0.08 ±0.08 


-0.07 ± 0.08 


lOg Sjr 


-0.09 ± 0.06 


-0.25 ± 0.07 


-0.24 ±0.07 


-0.32 ±0.07 


l.±0 


-0.03 ± 0.07 


0.37 ±0.07 


-0.07 ± 0.08 


0.23 ±0.07 


0.03 ±0.06 


logTi llg 


-0.01 ±0.06 


0.09 ±0.07 


0.15 ±0.07 


-0.18 ±0.08 


-0.03 ±0.07 


l.±0 


0.24 ±0.08 


-0.04 ±0.06 


-0.28 ±0.06 


0.33 ±0.09 


log 1ZE pk 


-0.15 ±0.06 


-0.21 ±0.06 


-0.25 ±0.06 


-0.37 ± 0.05 


0.37 ± 0.07 


0.24 ±0.08 


l.±0 


-0.03 ±0.07 


0.04 ±0.07 


-0.01 ±0.08 


log^ 


0.50 ± 0.05 


0.61 ±0.04 


0.61 ±0.04 


0.33 ±0.06 


-0.07 ± 0.08 


-0.04 ± 0.06 


-0.03 ±0.07 


l.±0 


0.58 ±0.04 


-0.20 ± 0.07 


log£pk 


0.24 ± 0.06 


0.14 ±0.08 


0.14 ±0.08 


0.08 ±0.08 


0.23 ±0.07 


-0.28 ± 0.06 


0.04 ±0.07 


0.58 ±0.04 


l.±0 


-0.28 ± 0.07 


a 


-0.26 ± 0.07 


-0.16 ±0.07 


-0.12 ±0.07 


-0.07 ± 0.08 


0.03 ± 0.06 


0.33 ±0.09 


-0.01 ±0.08 


-0.20 ±0.07 


-0.28 ±0.07 


l.±0 



TABLE 3 
Principal Components 



PCs log T 90 log T 50 log t log V log Sjr log Tlag log TZE pk log T log £ pk a Var % CVar % 



PC] 


-0.41 ± 0.02 


-0.49 ±0.01 


-0.45 ± 0.02 


-0.31 ± 0.03 


0.17 ±0.05 


0.04 ± 0.05 


0.18 ± 0.04 


-0.41 ±0.02 


-0.20 ± 0.05 


0.16 ± 0.05 


35.6 ± 1.6 


35.6 ± 1.6 


PC 2 


0.09 ± 0.08 


-0.11 ±0.08 


-0.13 ± 0.09 


-0.22 ±0.11 


0.48 ± 0.08 


-0.27 ± 0.25 


0.34 ±0.16 


0.25 ± 0.06 


0.56 ±0.09 


-0.35 ±0.14 


17. 4± 1.6 


53.0 ±5.1 


p c 3 


-0.17 ±0.07 


-0.18 ±0.05 


-0.21 ± 0.05 


0.29 ±0.11 


-0.22 ±0.19 


-0.66 ±0.10 


-0.46 ±0.17 


-0.13 ±0.11 


0.09 ± 0.20 


-0.29 ±0.15 


14.9 ± 1.4 


67.8 ± 5.5 


PC 4 


-0.30 ±0.14 


-0.09 ± 0.07 


-0.01 ± 0.06 


0.33 ± 0.16 


0.23 ± 0.20 


-0.02 ± 0.11 


-0.16 ±0.17 


0.28 ±0.10 


0.37 ±0.12 


0.71 ± 0.12 


8.7 ±0.8 


76.6 ± 2.1 


p c 5 


0.43 ± 0.13 


0.12 ±0.08 


0.06 ± 0.07 


-0.27 ± 0.25 


0.55 ±0.19 


-0.22 ±0.13 


-0.46 ±0.17 


-0.28 ±0.11 


-0.16 ± 0.17 


0.23 ± 0.23 


6.9 ± 0.8 


83.5 ± 3.1 


PCe 


-0.10 ± 0.19 


0.18 ±0.07 


0.04 ±0.10 


0.64 ± 0.14 


0.45 ± 0.22 


-0.14 ± 0.16 


0.38 ± 0.20 


-0.18 ±0.13 


-0.36 ± 0.12 


-0.12 ±0.20 


6.1 ±0.8 


89.6 ± 2.9 


PC 7 


-0.32 ±0.16 


-0.12 ±0.06 


0.16±0.17 


0.12 ± 0.19 


0.33 ±0.16 


0.55 ±0.13 


-0.46 ±0.13 


-0.07 ±0.13 


0.14±0.15 


-0.44 ±0.10 


4.5 ±0.5 


94.0 ± 1.1 


PC 8 


0.52 ± 0.1 1 


-0.03 ± 0.06 


-0.58 ± 0.09 


0.36 ± 0.10 


-0.09 ±0.13 


0.33 ±0.16 


-0.01 ±0.16 


-0.25 ±0.13 


0.29 ±0.11 


0.02 ± 0.13 


3.4 ±0.3 


97.4 ±0.6 


PC 9 


0.11 ± 0.10 


-0.15 ±0.05 


-0.42 ±0.10 


0.02 ± 0.08 


0.13 ±0.06 


0.06 ± 0.09 


-0.17 ±0.07 


0.70 ± 0.05 


-0.49 ± 0.06 


-0.08 ± 0.06 


2.1 ±0.2 


99.5 ±0.1 


^10 


0.36 ± 0.03 


-0.79 ±0.01 


0.45 ± 0.03 


0.19 ± 0.03 


-0.02 ± 0.02 


-0.01 ± 0.02 


0.09 ± 0.02 


0.01 ±0.03 


-0.06 ± 0.03 


0.02 ± 0.02 


0.5 ± 0.1 


100. ± 0. 



parameter, while the value of 7^' is derived from r' and the 
number of pulses. 
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